#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <cuda.h>

#define NUM_SMS (400)
#define NUM_THREADS_PER_SM (1536)
#define NUM_THREADS_PER_BLOCK (1024)
#define NUM_BLOCKS ((NUM_THREADS_PER_SM / NUM_THREADS_PER_BLOCK) * NUM_SMS)
#define NUM_ITERATIONS 5

// 128 MAD instructions
#define FMAD128(a, b) \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \
      a = b * a + b; \
      b = a * b + a; \

__shared__ float result[NUM_THREADS_PER_BLOCK];

__global__ void gflops()
{
    float a = result[threadIdx.x];  // this ensures the mads don't get compiled out
    float b = 1.01f;

    for (int i = 0; i < NUM_ITERATIONS; i++)
    {
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
        FMAD128(a, b);
    }
    result[threadIdx.x] = a + b;
}

int
main(int argc, char** argv) 
{
    // warmup
    gflops<<<NUM_BLOCKS, NUM_THREADS_PER_BLOCK>>>();
    cudaThreadSynchronize();

    cudaEvent_t start, stop; 
	float time; 
	cudaEventCreate(&start); 
	cudaEventCreate(&stop); 
	cudaEventRecord( start, 0 );

    gflops<<<NUM_BLOCKS, NUM_THREADS_PER_BLOCK>>>();
    cudaThreadSynchronize();
	cudaEventRecord( stop, 0 ); 
	cudaEventSynchronize( stop ); 
	cudaEventElapsedTime( &time, start, stop ); 
	cudaEventDestroy( start ); 
	cudaEventDestroy( stop );

    // output results
    printf( "Time: %f (ms)\n", time);
   // const double flops = 128 * 2 * 16 * NUM_ITERATIONS * (NUM_BLOCKS * NUM_THREADS_PER_BLOCK);
    
    //printf("Gflops: %f\n", (flops / (time / 1000.0f)) / 1e9 );
    const double flops = ((128 * 2 * 16)) * NUM_ITERATIONS * (NUM_BLOCKS * NUM_THREADS_PER_BLOCK / 1e9);
    
    printf("Gflops: %f\n", (flops / (time / 1000.0f) ));


	return 0;
}